
######## Set working directory to your local folder ######### 
setwd("C:/Users/brett/Google Drive/PFD and Crime/Paper/REStat/SubmissionFiles/Data and Code")



######## Custom functions ######################

copy_to_clipboard = function(x,sep="\t",col.names=T,...) { 
  write.table(x
              ,file = pipe("pbcopy")
              ,sep=sep
              ,col.names = col.names
              ,row.names = F
              ,quote = F,...)
}

nextweekday <- function(date, wday) {
  date <- as.Date(date)
  diff <- wday - wday(date)
  if( diff < 0 )
    diff <- diff + 7
  return(date + diff)
}

flattenlist <- function(x){  
  morelists <- sapply(x, function(xprime) class(xprime)[1]=="list")
  out <- c(x[!morelists], unlist(x[morelists], recursive=FALSE))
  if(sum(morelists)){ 
    Recall(out)
  }else{
    return(out)
  }
}

testFunction <- function (date_in) {
  return(tryCatch(as.Date(date_in), error=function(e) NULL))
}


theme_Publication <- function(base_size=14, base_family="helvetica") {
  library(grid)
  library(ggthemes)
  (theme_foundation(base_size=base_size, base_family=base_family)
    + theme(plot.title = element_text(face = "bold",
                                      size = rel(1.2), hjust = 0.5),
            text = element_text(),
            panel.background = element_rect(colour = NA),
            plot.background = element_rect(colour = NA),
            panel.border = element_rect(colour = NA),
            axis.title = element_text(face = "bold",size = rel(1)),
            axis.title.y = element_text(angle=90,vjust =2),
            axis.title.x = element_text(vjust = -0.2),
            axis.text = element_text(), 
            axis.line = element_line(colour="black"),
            axis.ticks = element_line(),
            panel.grid.major = element_line(colour="#f0f0f0"),
            panel.grid.minor = element_blank(),
            legend.key = element_rect(colour = NA),
            legend.position = "bottom",
            legend.direction = "horizontal",
            legend.key.size= unit(0.2, "cm"),
            legend.margin = unit(0, "cm"),
            legend.title = element_text(face="italic"),
            plot.margin=unit(c(10,5,5,5),"mm"),
            strip.background=element_rect(colour="#f0f0f0",fill="#f0f0f0"),
            strip.text = element_text(face="bold")
    ))
  
}

scale_fill_Publication <- function(...){
  library(scales)
  discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462","#7fc97f","#ef3b2c","#662506","#a6cee3","#fb9a99","#984ea3","#ffff33")), ...)
  
}

scale_colour_Publication <- function(...){
  library(scales)
  discrete_scale("colour","Publication",manual_pal(values = c("#386cb0","#fdb462","#7fc97f","#ef3b2c","#662506","#a6cee3","#fb9a99","#984ea3","#ffff33")), ...)
  
}

# Load Required Packages

library(viridis)
library(circular)
library(lubridate)
library(dplyr)
library(ggplot2)
library(chron)
library(timeDate)
library(readxl)
library(reshape2)
library(corrgram)
library(glmnet)
library( RcppBDT)
library(timeDate)
library(tis)
library(ggrepel)
library(stargazer)
library(broom)
library(scales)
library(multiwayvcov)
library(lmtest)
library(reporttools)
require(gridExtra)
library(reshape2)
library(formula.tools)
library(functional)
library(sandwich)
library(readr)
library(interplot)
library(margins)
library(scales)
library(car)
library(multcomp)
library(msm)









######## Loading Data ######### 

# Clean up

temp = paste(getwd(),"/APD Crime/",list.files(path="./APD Crime",pattern="*.csv"),sep="")
crime=do.call(rbind, lapply(temp, read.csv, header=T,stringsAsFactors=F))
crime$Received_Date=as.Date(crime$Received_Date,"%m/%d/%Y")
crime$Received_Time=unlist(lapply(crime$Received_Time,function(x) paste(paste(rep(0,6-nchar(x)),collapse = ""),x,sep="")))


# Load APD catagorization of event types:
#https://www.muni.org/Departments/police/Documents/PERF%20APD%20Deployment%20Study.pdf
MajorCADType <- read.csv("./MajorCADType.csv")[,1:2]
names(MajorCADType)=c("CallCode","MjrCADCat")

#Load Weather Data from NOAA for the airport
TedStevsAirportStationDailyWeather_NOAA_NCEI <- read.csv("./TedStevsAirportStationDailyWeather_NOAA_NCEI.csv")
TedStevsAirportStationDailyWeather_NOAA_NCEI$DATE=as.Date(TedStevsAirportStationDailyWeather_NOAA_NCEI$DATE,"%m/%d/%Y")


#Load Snap Payments
SnapPayments=read.csv("./SnapPaymentAgg.csv")[,-1]
SnapPaymentsD=read.csv("./SnapPaymentAgg.csv")
SnapPaymentsD[order(SnapPaymentsD$Persons),]

Inflation <- read.csv("./Inflation.csv")
Inflation=melt(Inflation,id.vars = "Year",variable.name = "Mon",value.name = "Inflation")
Inflation$cpi=Inflation$Inflation/241.432
SnapPayments=merge(SnapPayments,Inflation,by=c("Year","Mon")) # about 10% of the Alaska population receives food stamps
SnapPayments$CostDec2016=SnapPayments$Costs/SnapPayments$cpi
SnapPayments=SnapPayments[,c("Year","Mon","CostDec2016","Persons","Costs")]
names(SnapPayments)[2]="Month"

# Load PFD dates and ammounts 
pfdDatesAmmounts <- read.csv("./pfdDatesAmmounts.csv")
pfdDatesAmmounts$First_DD_Date=as.Date(pfdDatesAmmounts$First_DD_Date,"%m/%d/%Y")
pfdDatesAmounts_checks <- read.csv("./pfdDatesAmounts_checks.csv")
pfdDatesAmounts_checks$PFDSend_Date=as.Date(pfdDatesAmounts_checks$First_DD_Date,"%m/%d/%Y")


#Clean Dates and Create some indicaters
crime$DayofYear=format(crime$Received_Date,"%m-%d")
crime$Month=format(crime$Received_Date,"%m")
crime$DayOfWeek=weekdays(as.Date(crime$Received_Date))


# Combine PFD Ammounts and Crime data
crime=merge(crime,pfdDatesAmmounts,by="Year")
crime$DaystoPFD=crime$First_DD_Date-crime$Received_Date


# Go from long events to wide counts
crimeCounts=summarize(group_by(crime,Received_Date,DaystoPFD,DayofYear,DayOfWeek,Month,Year,`Call.Type.Desc`),CrimeCount=length(Call_No))
crimeCountsWide=dcast(crimeCounts,Received_Date+DayofYear+DayOfWeek+DaystoPFD+Month+Year~`Call.Type.Desc`,fill=0)
crimeCountsWide=merge(crimeCountsWide,pfdDatesAmmounts,by="Year")

crimeCountsWide$AfterPFD=0
crimeCountsWide$AfterPFD[which(crimeCountsWide$DaystoPFD<=0)]=1
crimeCountsWide$AfterPFD=factor(crimeCountsWide$AfterPFD)

# Add APD call codes (CAD) by major activity type (call for service CFS, Self-initiated activity SIA, or Admin)
crime_forCADmrg=crime[,c("Call.Type.Desc","Received_Date","DaystoPFD","Call_No")]
crime_forCADmrg$Call.Type.Desc=make.names(crime_forCADmrg$Call.Type.Desc)
crime_forCADmrg=merge(crime_forCADmrg,MajorCADType,by.x="Call.Type.Desc",by.y="CallCode")
crime_forCADmrg=summarize(group_by(crime_forCADmrg, Received_Date, DaystoPFD, MjrCADCat), CrimeCount=length(Call_No))
crime_forCADmrg_wide=dcast(crime_forCADmrg,Received_Date+DaystoPFD~MjrCADCat)
names(crime_forCADmrg_wide)[c(-1,-2)]=paste("MjrCADCd_",names(crime_forCADmrg_wide)[c(-1,-2)],sep="")

crimeCountsWide=merge(crimeCountsWide,crime_forCADmrg_wide[,c("Received_Date","MjrCADCd_Admin", "MjrCADCd_CFS", "MjrCADCd_SIA")], all.x=T)


### Find first weekends every month
crimeCountsWide2=crimeCountsWide
crimeCountsWide2$weekofyear=week(crimeCountsWide2$Received_Date)
crimeCountsWide2$dayofmonth=day(crimeCountsWide2$Received_Date)
crimeCountsWide2$dayofyear=yday(crimeCountsWide2$Received_Date)

crimeCountsWide2=mutate(group_by(crimeCountsWide2,DayOfWeek,Month,Year),FirstDayofMonthForDayOfWeek=min(dayofmonth))
crimeCountsWide2$FirstDOWmonth="NotFirstXofMonth"
crimeCountsWide2[which(crimeCountsWide2$FirstDayofMonthForDayOfWeek==crimeCountsWide2$dayofmonth),"FirstDOWmonth"]="FirstXofMonth"

### Day of week factor levels
crimeCountsWide2$DayOfWeek=factor(crimeCountsWide2$DayOfWeek,levels=c("Sunday","Monday", "Tuesday","Wednesday","Thursday","Friday","Saturday"))
crimeCountsWide2$DayOfWeekShort=substring(crimeCountsWide2$DayOfWeek,1,3)
crimeCountsWide2$DayOfWeekShort=factor(crimeCountsWide2$DayOfWeekShort,levels=c("Sun","Mon", "Tue","Wed","Thu","Fri","Sat"))

### Branch event data for enforcement effect analysis
crimetd=crime
crimetd$Date_Time=as.POSIXct(strptime(paste(crimetd$Received_Date,crimetd$Received_Time),"%Y-%m-%d %H%M%S"))
crimetd$Received_Hour=hour(crimetd$Date_Time)

crimetd$Call.Type.Desc=make.names(crimetd$Call.Type.Desc)
crimetd=merge(crimetd,MajorCADType,by.x="Call.Type.Desc",by.y="CallCode")


crimetd_sb=subset(crimetd,MjrCADCat!="")
crimetd_sb=summarise(group_by(crimetd_sb,Received_Date,Received_Hour,DaystoPFD,Year,Month,DayOfWeek,MjrCADCat),Incidents=length(Call_No))
crimetd_sb=dcast(crimetd_sb,Received_Date+Received_Hour+DaystoPFD+Year+Month+DayOfWeek~MjrCADCat,fill=0)
crimetd_sb$SIA_CFS_ratio=crimetd_sb$SIA/crimetd_sb$CFS
crimetd_sb$SIA_CFS_ratio[!is.finite(crimetd_sb$SIA_CFS_ratio)]=NA
crimetd_sb$DoubleStaffed=crimetd_sb$Received_Hour%in%c(0,7,8,15,16,23)

crimetd_sb=crimetd_sb[order(crimetd_sb$Received_Date,crimetd_sb$Received_Hour),]
crimetd_sb$SIA_CFS_ratioL1=dplyr::lag(crimetd_sb$SIA_CFS_ratio)
crimetd_sb$lnSIA=log(crimetd_sb$SIA)
crimetd_sb$lnSIA[!is.finite(crimetd_sb$lnSIA)]=NA
crimetd_sb$lnSIA_L1=dplyr::lag(crimetd_sb$lnSIA)

crimetd_sb$dayofmonth=day(crimetd_sb$Received_Date)
crimetd_sb$IsFederalHoliday=is.element(format(crimetd_sb$Received_Date,"%Y%m%d"),federalHolidays(2000:2016,board=T))
crimetd_sb=merge(crimetd_sb,TedStevsAirportStationDailyWeather_NOAA_NCEI,by.x="Received_Date",by.y="DATE",all.x=T)



######## Data Manipulation, Create Crime wide 3 ######### 

crimeCountsWide3= crimeCountsWide2
names(crimeCountsWide3)=make.names(names(crimeCountsWide3))
crimeCountsWide3$DateOrder=as.numeric(crimeCountsWide3$Received_Date)-10957
crimeCountsWide3$MonthID=as.numeric(crimeCountsWide3$Month)


#Define holidays - Floating
crimeCountsWide3$IsFederalHoliday=is.element(format(crimeCountsWide3$Received_Date,"%Y%m%d"),federalHolidays(2000:2016,board=T))
crimeCountsWide3$HalloweenOrWkndBF=is.element(format(crimeCountsWide3$Received_Date,"%Y%m%d"),c(paste(2000:2016,"1031",sep=""), as.character(latestSaturday(as.Date(paste(2000:2016,"10-31",sep="-")))),as.character(latestFriday(as.Date(paste(2000:2016,"10-31",sep="-")))) ) )
crimeCountsWide3$ColumbusOrWnknd=is.element(format(crimeCountsWide3$Received_Date,"%Y%m%d"),c(federalHolidays(2000:2016,board=T)[which(names(federalHolidays(2000:2016,board=T))=="Columbus")],federalHolidays(2000:2016,board=T)[which(names(federalHolidays(2000:2016,board=T))=="Columbus")]-1,federalHolidays(2000:2016,board=T)[which(names(federalHolidays(2000:2016,board=T))=="Columbus")]-2))
crimeCountsWide3$LaborOrWnknd=is.element(format(crimeCountsWide3$Received_Date,"%Y%m%d"),c(federalHolidays(2000:2016,board=T)[which(names(federalHolidays(2000:2016,board=T))=="Labor")],federalHolidays(2000:2016,board=T)[which(names(federalHolidays(2000:2016,board=T))=="Labor")]-1,federalHolidays(2000:2016,board=T)[which(names(federalHolidays(2000:2016,board=T))=="Labor")]-2))
crimeCountsWide3$Thanksgiving=is.element(format(crimeCountsWide3$Received_Date,"%Y%m%d"),c(federalHolidays(2000:2016,board=T)[which(names(federalHolidays(2000:2016,board=T))=="Thanksgiving")]))



#Define holidays - Rigid
crimeCountsWide3$JulyForth=crimeCountsWide3$DayofYear=="07-04"
crimeCountsWide3$Christmas=crimeCountsWide3$DayofYear=="12-25"
crimeCountsWide3$NewYears=crimeCountsWide3$DayofYear=="12-31" | crimeCountsWide3$DayofYear=="01-01"
crimeCountsWide3$StPat=crimeCountsWide3$DayofYear=="03-17"
crimeCountsWide3$CincoDeMayo=crimeCountsWide3$DayofYear=="05-04"
crimeCountsWide3$Iditarod =crimeCountsWide3$FirstDOWmonth=="FirstXofMonth" & crimeCountsWide3$DayOfWeekShort=="Sat" & crimeCountsWide3$Month=="03"
crimeCountsWide3$SuperBowl=is.element(format(crimeCountsWide3$Received_Date,"%Y-%m-%d"),c("1999-01-31","2000-01-30","2001-01-28","2002-02-03","2003-01-26","2004-02-01","2005-02-06","2006-02-05","2007-02-04","2008-02-03","2009-02-01","2010-02-07","2011-02-06","2012-02-05","2013-02-03","2014-02-02","2015-02-01","2016-02-07","2017-02-05"))


# Other Payouts: Military Pay day
crimeCountsWide3$MilitaryPayDay=0
crimeCountsWide3$MilitaryPayDay_dayAfter=0
crimeCountsWide3[which(is.element(crimeCountsWide3$Received_Date, c(crimeCountsWide3$Received_Date[which(is.element(crimeCountsWide3$dayofmonth,c(1,15))& isBusinessDay(crimeCountsWide3$Received_Date))], as.Date(previousBusinessDay(crimeCountsWide3$Received_Date[which(is.element(crimeCountsWide3$dayofmonth,c(1,15))& !isBusinessDay(crimeCountsWide3$Received_Date))]))))),"MilitaryPayDay"]=1

crimeCountsWide3[which(is.element(crimeCountsWide3$Received_Date, c(crimeCountsWide3$Received_Date[which(is.element(crimeCountsWide3$dayofmonth,c(1,15))& isBusinessDay(crimeCountsWide3$Received_Date))], as.Date(previousBusinessDay(crimeCountsWide3$Received_Date[which(is.element(crimeCountsWide3$dayofmonth,c(1,15))& !isBusinessDay(crimeCountsWide3$Received_Date))])) )+1 )),"MilitaryPayDay_dayAfter"]=1


crimeCountsWide3$MilitaryPayDay_ORdayAfter=0
crimeCountsWide3$MilitaryPayDay_ORdayAfter=crimeCountsWide3$MilitaryPayDay+crimeCountsWide3$MilitaryPayDay_dayAfter


# Other Payouts: Social Security
crimeCountsWide3$SSwed=0
tS = timeSequence(from = "2000-01-01", to = "2016-12-31", by = "month")
crimeCountsWide3[which(is.element(crimeCountsWide3$Received_Date,as.Date(timeNthNdayInMonth(tS, nday = 3, nth = 2)))),"SSwed"]=1
crimeCountsWide3[which(is.element(crimeCountsWide3$Received_Date,as.Date(timeNthNdayInMonth(tS, nday = 3, nth = 3)))),"SSwed"]=1
crimeCountsWide3[which(is.element(crimeCountsWide3$Received_Date,as.Date(timeNthNdayInMonth(tS, nday = 3, nth = 4)))),"SSwed"]=1


# Relative to PFD time dummies
crimeCountsWide3 $PFDweekend=0
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==0),"PFDweekend"]=1
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-1),"PFDweekend"]=1
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-2),"PFDweekend"]=1
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-3),"PFDweekend"]=1

crimeCountsWide3$TimeAfterPFDDummy="1NoPFD"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==0),"TimeAfterPFDDummy"]="2FirstWeekend"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-1),"TimeAfterPFDDummy"]="2FirstWeekend"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-2),"TimeAfterPFDDummy"]="2FirstWeekend"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-3),"TimeAfterPFDDummy"]="2FirstWeekend"

crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-4),"TimeAfterPFDDummy"]="3FirstWeekAfter"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-5),"TimeAfterPFDDummy"]="3FirstWeekAfter"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-6),"TimeAfterPFDDummy"]="3FirstWeekAfter"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-7),"TimeAfterPFDDummy"]="3FirstWeekAfter"

crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-8),"TimeAfterPFDDummy"]="4SecondWeekendAfter"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-9),"TimeAfterPFDDummy"]="4SecondWeekendAfter"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-10),"TimeAfterPFDDummy"]="4SecondWeekendAfter"

crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-11),"TimeAfterPFDDummy"]="5SecondWeekAfter"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-12),"TimeAfterPFDDummy"]="5SecondWeekAfter"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-13),"TimeAfterPFDDummy"]="5SecondWeekAfter"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-14),"TimeAfterPFDDummy"]="5SecondWeekAfter"

crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-15),"TimeAfterPFDDummy"]="6ThirdWeekendAfter"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-16),"TimeAfterPFDDummy"]="6ThirdWeekendAfter"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-17),"TimeAfterPFDDummy"]="6ThirdWeekendAfter"

crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-18),"TimeAfterPFDDummy"]="7ThirdWeekAfter"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-19),"TimeAfterPFDDummy"]="7ThirdWeekAfter"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-20),"TimeAfterPFDDummy"]="7ThirdWeekAfter"
crimeCountsWide3[which(crimeCountsWide3 $DaystoPFD ==-21),"TimeAfterPFDDummy"]="7ThirdWeekAfter"

# Add weather data
crimeCountsWide3=merge(crimeCountsWide3,TedStevsAirportStationDailyWeather_NOAA_NCEI,by.x="Received_Date",by.y="DATE")

# "Fixed effects" label for reference
crimeCountsWide3$MonthFixedEffect=crimeCountsWide3$Month
crimeCountsWide3$YearFixedEffect=crimeCountsWide3$Year

#Year-week for clustering
crimeCountsWide3$YearWeek=paste(crimeCountsWide3$Year,crimeCountsWide3$weekofyear,sep="-")

# Windows around Treatment
crimeCountsWide3$WeekBeforeWindow=crimeCountsWide3$DaystoPFD>0 & crimeCountsWide3$DaystoPFD<=6
crimeCountsWide3$WeekAftereWindow=crimeCountsWide3$DaystoPFD<0 & crimeCountsWide3$DaystoPFD>=-7

# True PFD Weekend sun=1, Sat=7
crimeCountsWide3$PFDweekend=0
crimeCountsWide3[which(is.element(crimeCountsWide3$Received_Date, nextweekday(crimeCountsWide3[which(crimeCountsWide3$DaystoPFD==0),"Received_Date"], 6))),"PFDweekend"]=1 # Friday
crimeCountsWide3[which(is.element(crimeCountsWide3$Received_Date, nextweekday(crimeCountsWide3[which(crimeCountsWide3$DaystoPFD==0),"Received_Date"], 7))),"PFDweekend"]=1 # Sat
crimeCountsWide3[which(is.element(crimeCountsWide3$Received_Date, nextweekday(crimeCountsWide3[which(crimeCountsWide3$DaystoPFD==0),"Received_Date"], 1))),"PFDweekend"]=1 # Sun
head(subset(crimeCountsWide3,PFDweekend==1))
crimeCountsWide3$NumPFDweekend=as.numeric(crimeCountsWide3$PFDweekend)

#Numeric Day after
crimeCountsWide3$NumPFDdayAfter=as.numeric(crimeCountsWide3$DaystoPFD==-1)

# Treatment effect through sunday
crimeCountsWide3$ThroughSunday=0
crimeCountsWide3[which(is.element(crimeCountsWide3$Received_Date, nextweekday(crimeCountsWide3[which(crimeCountsWide3$DaystoPFD==0),"Received_Date"], 1))),"ThroughSunday"]=1 # Sun
crimeCountsWide3[which(is.element(crimeCountsWide3$Received_Date, nextweekday(crimeCountsWide3[which(crimeCountsWide3$DaystoPFD==0),"Received_Date"], 7))),"ThroughSunday"]=1 # Sat
crimeCountsWide3[which(is.element(crimeCountsWide3$Received_Date, nextweekday(crimeCountsWide3[which(crimeCountsWide3$DaystoPFD==0),"Received_Date"], 6)) & crimeCountsWide3$DaystoPFD<0 ),"ThroughSunday"]=1 # Friday
crimeCountsWide3[which(is.element(crimeCountsWide3$Received_Date, nextweekday(crimeCountsWide3[which(crimeCountsWide3$DaystoPFD==0),"Received_Date"], 5)) & crimeCountsWide3$DaystoPFD<0 ),"ThroughSunday"]=1 # Thursday
crimeCountsWide3[which(is.element(crimeCountsWide3$Received_Date, nextweekday(crimeCountsWide3[which(crimeCountsWide3$DaystoPFD==0),"Received_Date"], 4)) & crimeCountsWide3$DaystoPFD<0 ),"ThroughSunday"]=1 # Wed


# Normalized PFD payments
crimeCountsWide3$PFD_FirstDD_Amount_FromMin=(crimeCountsWide3$Amount_First_Date-min(crimeCountsWide3$Amount_First_Date))/100000000
crimeCountsWide3$PFD_FirstDD_Amount_FromAvg=(crimeCountsWide3$Amount_First_Date-mean(crimeCountsWide3$Amount_First_Date))/100000000
crimeCountsWide3$Amount_First_Date_100m=crimeCountsWide3$Amount_First_Date/100000000

crimeCountsWide3$PFD_Amount_FromMin=(crimeCountsWide3$PFD_Amount-min(crimeCountsWide3$PFD_Amount))
crimeCountsWide3$PFD_Amount_FromAvg=(crimeCountsWide3$PFD_Amount-mean(unique(crimeCountsWide3$PFD_Amount)))/100
crimeCountsWide3$PFD_Amount_100dol=crimeCountsWide3$PFD_Amount/1000

crimeCountsWide3=merge(crimeCountsWide3,summarize(group_by(crimeCountsWide3[crimeCountsWide3$HalloweenOrWkndBF,c("Year","DaystoPFD")],Year),DaysToHolloween=abs(max(DaystoPFD))),by="Year")
crimeCountsWide3$PFDtoHalloween=(crimeCountsWide3$DaystoPFD*-1<=crimeCountsWide3$DaysToHolloween) & crimeCountsWide3$DaystoPFD<=0

crimeCountsWide3$NumWeekAftereWindow=as.numeric(crimeCountsWide3$WeekAftereWindow)
crimeCountsWide3$NumPFDtoHalloween=as.numeric(crimeCountsWide3$PFDtoHalloween)
crimeCountsWide3$NumPFDtoHalloween=as.numeric(crimeCountsWide3$PFDtoHalloween)
crimeCountsWide3$Num2Weeks=as.numeric(crimeCountsWide3$DaystoPFD<0 & crimeCountsWide3$DaystoPFD>=-14)

#Conventional Persistance
crimeCountsWide3$NumDay8to14=as.numeric(crimeCountsWide3$DaystoPFD<(-8) & crimeCountsWide3$DaystoPFD>=-14)
crimeCountsWide3$NumDay15to21=as.numeric(crimeCountsWide3$DaystoPFD<(-14) & crimeCountsWide3$DaystoPFD>=-21)
crimeCountsWide3$NumDay22to28=as.numeric(crimeCountsWide3$DaystoPFD<(-22) & crimeCountsWide3$DaystoPFD>=-28)

crimeCountsWide3$PFD_Amount_Persons_FrstDay_100k=(crimeCountsWide3$Amount_First_Date/crimeCountsWide3$PFD_Amount)/100000

#add holiday label for easy reference
names(crimeCountsWide3)[which(names(crimeCountsWide3)%in%c("NewYears", "SuperBowl", "Iditarod", "StPat", "CincoDeMayo", "JulyForth", "LaborOrWnknd", "ColumbusOrWnknd", "HalloweenOrWkndBF", "Thanksgiving", "Christmas", "IsFederalHoliday") )]=paste("Holiday_",names(crimeCountsWide3)[which(names(crimeCountsWide3)%in%c("NewYears", "SuperBowl", "Iditarod", "StPat", "CincoDeMayo", "JulyForth", "LaborOrWnknd", "ColumbusOrWnknd", "HalloweenOrWkndBF", "Thanksgiving", "Christmas", "IsFederalHoliday") )],sep="")




########### Results ########### 


###########** Results: One day after ########### 

FormulaList=list()
FormulBase=as.formula(DRUNK.PROBLEM~PFDweekend)


# Day After
rhs(FormulBase)=quote(NumPFDdayAfter+MilitaryPayDay_ORdayAfter+Holiday_NewYears+Holiday_SuperBowl+Holiday_Iditarod+Holiday_StPat+Holiday_CincoDeMayo+Holiday_JulyForth+Holiday_LaborOrWnknd+Holiday_ColumbusOrWnknd+Holiday_HalloweenOrWkndBF+Holiday_Thanksgiving+Holiday_Christmas+Holiday_IsFederalHoliday+poly(PRCP,3)+poly(SNWD,3)+poly(TMAX,3)+DayOfWeek + poly(dayofmonth,5) +factor(MonthFixedEffect)*factor(YearFixedEffect))


# Name LHS
FormulaList[["Violence"]]=FormulBase
FormulaList[["SubstancePart"]]=FormulBase
FormulaList[["SubstanceFull"]]=FormulBase
FormulaList[["Property"]]=FormulBase
FormulaList[["Party"]]=FormulBase
FormulaList[["MedicalAssist"]]=FormulBase


DesiredOutcomes=c("Violence","SubstanceFull","Property","Party","MedicalAssist")
#Define LHS
lhs(FormulaList[["SubstancePart"]])=quote(DRIVING.WHILE.INTOXICATED+DRUGS.FORGED.PERSCRIPTION+DRUNK.TRANSPORT+DRUNK.PROBLEM)
lhs(FormulaList[["SubstanceFull"]])=quote(DRIVING.WHILE.INTOXICATED+HIT.AND.RUN+HIT.AND.RUN.W.INJURY+DRUGS.FORGED.PERSCRIPTION+DRUNK.TRANSPORT+DRUNK.PROBLEM+LIQUOR.LAW.VIOLATION)
lhs(FormulaList[["Party"]])=quote(LOUD.DISRUPTIVE.PARTY+NOISE.VIOLATION)
lhs(FormulaList[["Violence"]])=quote(SEXUAL.ASLT.IN.PROGRESS+SEXUAL.ASSAULT.OF.MINOR+SEXUAL.ASSAULT+ASSAULT+ASSAULT.WITH.A.WEAPON+HOMICIDE)
lhs(FormulaList[["Property"]])=quote(BURGLARY+THEFT+SHOPLIFTER+STOLEN.PROPERTY+STOLEN.VEHICLE+ROBBERY+STRONGARM.ROBBERY+BURGLARY.INPROGRESS)
lhs(FormulaList[["MedicalAssist"]])=quote(MEDIC.ASSIST)



FormulaList_MAIN=FormulaList

# Run models, store in list
ModelList=lapply(FormulaList,function(x) lm(x,subset(crimeCountsWide3, DaystoPFD!=0)))
# Apply cluster robust method
DayAfterModels_Basic=ModelList
DayAfterModelErrors_Basic_VCOV=lapply(ModelList, function(x) NeweyWest(x)) 
ClstrdErr=lapply(DayAfterModelErrors_Basic_VCOV, function(x) sqrt(diag(x)) )
DayAfterModelErrors_Basic=ClstrdErr



##### *****Table 2 ######

# Desired outcomes for main text table
DesiredOutcomes=c("Violence","SubstanceFull","Property","Party","MedicalAssist")

resultPvalues=unlist(mapply(function(x,vcov) tidy(coeftest(x,vcov=vcov))[which(tidy(x)$term=="NumPFDdayAfter"),"p.value"],x=DayAfterModels_Basic[DesiredOutcomes],vcov=DayAfterModelErrors_Basic_VCOV[DesiredOutcomes]))


stargazer(DayAfterModels_Basic[DesiredOutcomes],omit=c("YearFixedEffect","poly","Holiday_","DayOfWeek","dayofmonth","MonthFixedEffect","factor.MonthFixedEffect...:factor.YearFixedEffect."),omit.labels=c("Year Fixed Effects","Weather","Holiday Effects","Day of Week Fixed Effects","Day-of-Month 5th Order Spline","Month Fixed Effects","Month x Year Effects"), column.labels = names(FormulaList[DesiredOutcomes]),se=DayAfterModelErrors_Basic[DesiredOutcomes],no.space=T,keep.stat=c("n","adj.rsq"),notes="Newey-West Robust Errors.") 


print(xtable(rbind(
  data.frame(variable="P-val",t(paste("[",format(resultPvalues,digits=1),"]",sep=""))),
  data.frame(variable="Bonferroni P-val",t(paste("[",format(p.adjust(resultPvalues,"bonferroni"),digits=2),"]",sep="")))
)), include.rownames=FALSE)


print(xtable(rbind(
  data.frame(variable="Mean Daily Call Count",round(t(unlist(lapply(FormulaList[DesiredOutcomes], function(x) mean(rowSums(subset(crimeCountsWide3)[,lhs.vars(x),drop=FALSE])) )) ),2)),
  data.frame(variable="St. dev Call Count",round(t(unlist(lapply(FormulaList[DesiredOutcomes], function(x) sd(rowSums(subset(crimeCountsWide3)[,lhs.vars(x),drop=FALSE])) ))),2)))), include.rownames=FALSE)


##### *****Appendix: Table A:2 ######
# Appendix Table, Full results, holidays, adjusted p-values

DesiredOutcomes=c("Violence","SubstancePart","SubstanceFull","Property","Party","MedicalAssist")

stargazer(DayAfterModels_Basic[DesiredOutcomes],omit=c("YearFixedEffect","poly","DayOfWeek","dayofmonth","MonthFixedEffect","factor.MonthFixedEffect...:factor.YearFixedEffect."),omit.labels=c("Year Fixed Effects","Weather","Day of Week Fixed Effects","Day-of-Month 5th Order Spline","Month Fixed Effects","Month x Year Effects"),covariate.labels=c(NA,"Mil. Pay Day/Day After","New Years Day/Eve", "Super Bowl", "Iditarod", "St Patricks Day", "Cinco de Mayo", "July 4th", "Labor Day Weekend", "Columbus Day Weekend", "Halloween and Weekend", "Thanksgiving", "Christmas", "Federal Holiday"), column.labels = names(FormulaList[DesiredOutcomes]),se=DayAfterModelErrors_Basic[DesiredOutcomes],no.space=T,keep.stat=c("n","adj.rsq"),notes="Newey-West Robust Errors.") 

resultPvalues=unlist(mapply(function(x,vcov) tidy(coeftest(x,vcov=vcov))[which(tidy(x)$term=="NumPFDdayAfter"),"p.value"],x=DayAfterModels_Basic[DesiredOutcomes],vcov=DayAfterModelErrors_Basic_VCOV[DesiredOutcomes]))

# Correction of p-values for multiple hypothesis testing - loop through methods:
sapply(p.adjust.methods[p.adjust.methods != "fdr"], FUN=function(meth) p.adjust(resultPvalues,method=meth))

# Print p value corrections:
print(xtable(rbind(
  data.frame(variable="Unadjusted P-value",t(paste("[",format(resultPvalues,digits=1),"]",sep=""))),
  data.frame(variable="Bonferroni P-value",t(paste("[",format(p.adjust(resultPvalues,"bonferroni"),digits=2),"]",sep="")))
)), include.rownames=FALSE)



########### **Results: One day after: parsimonious controls ########### 


FormulaList=list()
FormulBase=as.formula(DRUNK.PROBLEM~PFDweekend)


# Day After
rhs(FormulBase)=quote(NumPFDdayAfter+DayOfWeek+factor(MonthFixedEffect)*factor(YearFixedEffect))

# Name LHS
FormulaList[["Violence"]]=FormulBase
FormulaList[["SubstancePart"]]=FormulBase
FormulaList[["SubstanceFull"]]=FormulBase
FormulaList[["Property"]]=FormulBase
FormulaList[["Party"]]=FormulBase
FormulaList[["MedicalAssist"]]=FormulBase


#Define LHS variables for the outcome variables of interest
lhs(FormulaList[["SubstancePart"]])=quote(DRIVING.WHILE.INTOXICATED+DRUGS.FORGED.PERSCRIPTION+DRUNK.TRANSPORT+DRUNK.PROBLEM)
lhs(FormulaList[["SubstanceFull"]])=quote(DRIVING.WHILE.INTOXICATED+HIT.AND.RUN+HIT.AND.RUN.W.INJURY+DRUGS.FORGED.PERSCRIPTION+DRUNK.TRANSPORT+DRUNK.PROBLEM+LIQUOR.LAW.VIOLATION)
lhs(FormulaList[["Party"]])=quote(LOUD.DISRUPTIVE.PARTY+NOISE.VIOLATION)
lhs(FormulaList[["Violence"]])=quote(SEXUAL.ASLT.IN.PROGRESS+SEXUAL.ASSAULT.OF.MINOR+SEXUAL.ASSAULT+ASSAULT+ASSAULT.WITH.A.WEAPON+HOMICIDE)
lhs(FormulaList[["Property"]])=quote(BURGLARY+THEFT+SHOPLIFTER+STOLEN.PROPERTY+STOLEN.VEHICLE+ROBBERY+STRONGARM.ROBBERY+BURGLARY.INPROGRESS)
lhs(FormulaList[["MedicalAssist"]])=quote(MEDIC.ASSIST)


# Run models, store in list
ModelList=lapply(FormulaList,function(x) lm(x,subset(crimeCountsWide3, DaystoPFD!=0)))
DayAfterModelVCOV_Parsimonious=lapply(ModelList, function(x) NeweyWest(x))
# Apply cluster robust method
ClstrdErr=lapply(DayAfterModelVCOV_Parsimonious, function(x) sqrt(diag(x)) )

DayAfterModels_Parsimonious=ModelList
DayAfterModelErrors_Parsimonious=ClstrdErr

# Desired outcomes for appendix table
DesiredOutcomes=c("Violence","SubstancePart","SubstanceFull","Property","Party","MedicalAssist")

##### *****Appendix: Table A:3 ######
stargazer(DayAfterModels_Parsimonious[DesiredOutcomes],omit=c("YearFixedEffect","poly","Holiday_","dayofmonth","MonthFixedEffect","factor.MonthFixedEffect...:factor.YearFixedEffect."),omit.labels=c("Year Fixed Effects","Weather","Holiday Effects","Day-of-Month 5th Order Spline","Month Fixed Effects","Month x Year Effects"), column.labels = names(FormulaList[DesiredOutcomes]),se=DayAfterModelErrors_Parsimonious[DesiredOutcomes],no.space=T,keep.stat=c("n","adj.rsq"),notes="Newey-West Robust Errors.") 


resultPvalues=unlist(mapply(function(x,vcov) tidy(coeftest(x,vcov=vcov))[which(tidy(x)$term=="NumPFDdayAfter"),"p.value"],x=DayAfterModels_Parsimonious[DesiredOutcomes],vcov=DayAfterModelVCOV_Parsimonious[DesiredOutcomes]))

# Print p value corrections:
print(xtable(rbind(
  data.frame(variable="Unadjusted P-value",t(paste("[",sprintf("%.3f",resultPvalues),"]",sep=""))),
  data.frame(variable="Bonferroni P-value",t(paste("[",sprintf("%.3f",p.adjust(resultPvalues,"bonferroni")),"]",sep="")))
)), include.rownames=FALSE)





###########** Results: One day after: Poisson ########### 


ModelList=lapply(FormulaList,function(x) glm(x,family="poisson",subset(crimeCountsWide3, DaystoPFD!=0)))
DayAfterModelVCOV_Poisson=lapply(ModelList, function(x) NeweyWest(x))
ClstrdErr=lapply(DayAfterModelVCOV_Poisson, function(x) sqrt(diag(x)) )
DayAfterModels_Poisson=ModelList
DayAfterModelErrors_Poisson=ClstrdErr


DesiredOutcomes=c("Violence","SubstancePart","SubstanceFull","Property","Party","MedicalAssist")

##### *****Appendix: Table A:4 ######
stargazer(DayAfterModels_Poisson[DesiredOutcomes],omit=c("YearFixedEffect","poly","DayOfWeek","dayofmonth","MonthFixedEffect","factor.MonthFixedEffect...:factor.YearFixedEffect."),omit.labels=c("Year Fixed Effects","Weather","Day of Week Fixed Effects","Day-of-Month 5th Order Spline","Month Fixed Effects","Month x Year Effects"),covariate.labels=c(NA,"Mil. Pay Day/Day After","New Years Day/Eve", "Super Bowl", "Iditarod", "St Patricks Day", "Cinco de Mayo", "July 4th", "Labor Day Weekend", "Columbus Day Weekend", "Halloween and Weekend", "Thanksgiving", "Christmas", "Federal Holiday"), column.labels = names(FormulaList[DesiredOutcomes]),se=DayAfterModelErrors_Poisson[DesiredOutcomes],no.space=T,keep.stat=c("n","adj.rsq"),notes="Newey-West Robust Errors.") 

resultPvalues=unlist(mapply(function(x,vcov) tidy(coeftest(x,vcov=vcov))[which(tidy(x)$term=="NumPFDdayAfter"),"p.value"],x=DayAfterModels_Poisson[DesiredOutcomes],vcov=DayAfterModelVCOV_Poisson[DesiredOutcomes]))

# Print p value corrections:
print(xtable(rbind(
  data.frame(variable="Unadjusted P-value",t(paste("[",format(resultPvalues,digits=1),"]",sep=""))),
  data.frame(variable="Bonferroni P-value",t(paste("[",format(p.adjust(resultPvalues,"bonferroni"),digits=2),"]",sep="")))
)), include.rownames=FALSE)

# Print mean and standard deviation of the outcome variables 
print(xtable(rbind(
  data.frame(variable="Mean Daily Call Count",round(t(unlist(lapply(FormulaList[DesiredOutcomes], function(x) mean(rowSums(subset(crimeCountsWide3)[,lhs.vars(x),drop=FALSE])) )) ),2)),
  data.frame(variable="St. dev Call Count",round(t(unlist(lapply(FormulaList[DesiredOutcomes], function(x) sd(rowSums(subset(crimeCountsWide3)[,lhs.vars(x),drop=FALSE])) ))),2)))), include.rownames=FALSE)






########### ** Disaggregate effects  ########### 

# Create vector of oucome variables
varlist_lhs=do.call(rbind,lapply(FormulaList,lhs.vars))
varlist_lhs_grp=unique(melt(varlist_lhs)[,c(1,3)])

# write.table(unique(melt(varlist_lhs)[,c(1,3)]), "clipboard", sep="\t", row.names=FALSE, col.names=FALSE)
Formula_Breakout=as.formula(DRUNK.PROBLEM~PFDweekend)
FormulaList_breakout=list()

breakout_LHS=unique(unlist(lapply(FormulaList,function(x) lhs.vars(x))))
breakout_LHS=make.names(unique(crime$Call.Type.Desc)[-106])
FormulaList_breakout=lapply(breakout_LHS,function(x) as.formula(paste(x,"~0")))


# Set RHS of breakout formula
rhs(FormulaList_breakout)=lapply(FormulaList_breakout,function(x) quote(NumWeekAftereWindow+MilitaryPayDay_ORdayAfter + Holiday_NewYears + Holiday_SuperBowl + 
                                                                          Holiday_Iditarod + Holiday_StPat + Holiday_CincoDeMayo + 
                                                                          Holiday_JulyForth + Holiday_LaborOrWnknd + Holiday_ColumbusOrWnknd + 
                                                                          Holiday_HalloweenOrWkndBF + Holiday_Thanksgiving + Holiday_Christmas + 
                                                                          Holiday_IsFederalHoliday + poly(PRCP, 3) + poly(SNWD, 3) + 
                                                                          poly(TMAX, 3) + DayOfWeek + poly(dayofmonth, 5) + factor(MonthFixedEffect) * 
                                                                          factor(YearFixedEffect)))

# Mean number of incidents for table
averageIncidents=colMeans(crimeCountsWide3[,breakout_LHS])
averageIncidents=data.frame(type=names(averageIncidents),averageIncidents)

# Run OLS model:
ModelList_breakout=lapply(FormulaList_breakout,function(x) data.frame(Outcome=as.character(lhs(x)),tidy(lm(x,subset(crimeCountsWide3,DaystoPFD!=0)))[2,]))
ModelList_breakoutClean=do.call(rbind,ModelList_breakout)
names(ModelList_breakoutClean)[-1]=paste("lm",names(ModelList_breakoutClean)[-1],sep=".")

# Run Poisson model:
ModelList_breakout_POISSON=lapply(FormulaList_breakout,function(x) return(tryCatch(data.frame(Outcome=as.character(lhs(x)),tidy(glm(x,subset(crimeCountsWide3,DaystoPFD!=0),family="poisson"))[2,]), error=function(e) NULL)) )
ModelList_breakoutClean_Pois=do.call(rbind,ModelList_breakout_POISSON)

# Merge OLS and Poisson
ModelList_breakoutClean_Pois=merge(ModelList_breakoutClean_Pois,averageIncidents,by.x="Outcome",by.y="type")
ModelList_breakoutClean_Pois=merge(ModelList_breakoutClean_Pois,varlist_lhs_grp,by.x="Outcome",by.y="value",all.x=T)
ModelList_breakoutClean_both=merge(ModelList_breakoutClean_Pois,ModelList_breakoutClean,by="Outcome")

# Multiple hypothesis test corrects:
ModelList_breakoutClean_both$adjustedP=NA
ModelList_breakoutClean_both$lm.adjustedP=NA
ModelList_breakoutClean_both[which(!is.na(ModelList_breakoutClean_both$Var1)),"adjustedP"]=p.adjust(ModelList_breakoutClean_both[which(!is.na(ModelList_breakoutClean_both$Var1)),"p.value"],method = "bonferroni")
ModelList_breakoutClean_both[which(!is.na(ModelList_breakoutClean_both$Var1)),"lm.adjustedP"]=p.adjust(ModelList_breakoutClean_both[which(!is.na(ModelList_breakoutClean_both$Var1)),"lm.p.value"],method = "bonferroni")

# Write out results table: 
write.csv(ModelList_breakoutClean_both,"./Results/ModelList_breakoutClean_both.csv")

#IncidentCountsByYear=summarise(group_by(crime,Year,Call.Type.Desc),AnnualCount=length(Call_No))
#IncidentCountsByYear$Call.Type.Desc=make.names(IncidentCountsByYear$Call.Type.Desc)
#write.csv(IncidentCountsByYear,"./Results/AnnualIncidentCounts.csv")


##### *****Appendix: Table A:5 ######
subset(ModelList_breakoutClean_both[order(ModelList_breakoutClean_both$Var1),],!is.na(Var1))



########### **Results: Check day after results in smaller window: 09-2016 ########### 

DesiredOutcomes=c("Violence","SubstanceFull","Property","Party","MedicalAssist")

# Run models, store in list
ModelList=lapply(FormulaList_MAIN[DesiredOutcomes],function(x) lm(x,subset(crimeCountsWide3,Year<=2009 & DaystoPFD!=0)))
# Apply cluster robust method
ClstrdErr=lapply(ModelList, function(x) sqrt(diag(NeweyWest(x))) )
DayAfterModels_Basic=ModelList
DayAfterModelErrors_Basic=ClstrdErr

Pre09Models=list()
Pre09Models[["Models"]]=ModelList
Pre09Models[["Errors"]]=ClstrdErr


# Run models, store in list
ModelList=lapply(FormulaList_MAIN[DesiredOutcomes],function(x) lm(x,subset(crimeCountsWide3,Year>2009 & DaystoPFD!=0)))
# Apply cluster robust method
ClstrdErr=lapply(ModelList, function(x) sqrt(diag(NeweyWest(x))) )
DayAfterModels_Basic=ModelList
DayAfterModelErrors_Basic=ClstrdErr
Post09Models=list()
Post09Models[["Models"]]=ModelList
Post09Models[["Errors"]]=ClstrdErr

##### *****Appendix: Table A:7 ######
stargazer(c(Pre09Models$Models,Post09Models$Models),omit=c("YearFixedEffect","poly","Holiday_","DayOfWeek","dayofmonth","MonthFixedEffect","factor.MonthFixedEffect...:factor.YearFixedEffect."),omit.labels=c("Year Fixed Effects","Weather","Holiday Effects","Day of Week Fixed Effects","Day-of-Month 5th Order Poly. Trend","Month Fixed Effects","Month x Year Effects"), column.labels = c(names(ModelList),names(ModelList)),se=c(Pre09Models$Errors,Post09Models$Errors),no.space=T,keep.stat=c("n","adj.rsq"),notes="Newey-West Robust Errors.") 







########### Appendix C: Enforcemnt effect:  ########### 


FormulaList_Enforcement=list()

FormulaList_Enforcement[["Admin"]]=FormulBase
FormulaList_Enforcement[["CFS"]]=FormulBase
FormulaList_Enforcement[["SIA"]]=FormulBase
FormulaList_Enforcement[["CFS/SIA"]]=FormulBase
FormulaList_Enforcement[["CFS/Admin"]]=FormulBase
FormulaList_Enforcement[["SIA/Admin"]]=FormulBase

lhs(FormulaList_Enforcement[["Admin"]])=quote(MjrCADCd_Admin)
lhs(FormulaList_Enforcement[["CFS"]])=quote(MjrCADCd_CFS)
lhs(FormulaList_Enforcement[["SIA"]])=quote(MjrCADCd_SIA)
lhs(FormulaList_Enforcement[["CFS/SIA"]])=quote(MjrCADCd_CFS/MjrCADCd_SIA)
lhs(FormulaList_Enforcement[["CFS/Admin"]])=quote(MjrCADCd_CFS/MjrCADCd_Admin)
lhs(FormulaList_Enforcement[["SIA/Admin"]])=quote(MjrCADCd_SIA/MjrCADCd_Admin)

FormulaList_Enforcement_hr=FormulaList_Enforcement[c("Admin","CFS","SIA")]
lhs(FormulaList_Enforcement_hr[["Admin"]])=quote(Admin)
lhs(FormulaList_Enforcement_hr[["CFS"]])=quote(CFS)
lhs(FormulaList_Enforcement_hr[["SIA"]])=quote(SIA)


rhs(FormulaList_Enforcement_hr)=lapply(FormulaList_Enforcement_hr,function(x) quote(I(DaystoPFD==-1)+poly(PRCP,3)+poly(SNWD,3)+poly(TMAX,3)+IsFederalHoliday+ DoubleStaffed+DayOfWeek+ poly(dayofmonth,5)+ poly(dayofmonth,5) +factor(Month)*factor(Year)))


ModelList_hourlyenfocement=lapply(FormulaList_Enforcement_hr,function(x) glm(x,family="poisson",subset(crimetd_sb, DaystoPFD!=0)))

#DayAfterModelVCOV_hourlyenfocement=lapply(ModelList_hourlyenfocement, function(x) NeweyWest(x))
#ClstrdErr_hourlyenfocement=lapply(DayAfterModelVCOV_hourlyenfocement, function(x) sqrt(diag(x)) )

unlist(lapply(FormulaList_Enforcement_hr,function(x) mean(crimetd_sb[,lhs.vars(x)])))


stargazer(ModelList_hourlyenfocement,omit=c("DayOfWeek","poly","dayofmonth","Month","Year"),omit.labels=c("Day of Week","Weather","Day of month 5th order poly. trend","","Month x Year FEs"), column.labels = names(FormulaList_Enforcement_hr),no.space=T,keep.stat=c("n","adj.rsq"),notes="Newey-West Robust Errors.") 

##### *****Appendix C: Table C:1 ######
stargazer(ModelList_hourlyenfocement,omit=c("DayOfWeek","poly","dayofmonth","Month","Year"),omit.labels=c("Day of Week","Weather","Day of month 5th order poly. trend","","Month x Year FEs"), column.labels = names(FormulaList_Enforcement_hr),no.space=T,keep.stat=c("n","adj.rsq"),notes="Newey-West Robust Errors.",type="text") 



# Daily counts
# Run models, store in list
ModelList_enforcement=lapply(FormulaList_Enforcement,function(x) lm(x,subset(crimeCountsWide3, DaystoPFD!=0)))
# Apply cluster robust method
DayAfterModelErrors_Enforce_VCOV=lapply(ModelList_enforcement, function(x) NeweyWest(x)) 
ClstrdErr_Enforce=lapply(DayAfterModelErrors_Enforce_VCOV, function(x) sqrt(diag(x)) )


stargazer(ModelList_enforcement,omit=c("YearFixedEffect","poly","DayOfWeek","dayofmonth","MonthFixedEffect","factor.MonthFixedEffect...:factor.YearFixedEffect.","weekofyear"),omit.labels=c("Year Fixed Effects","Weather","Day of Week Fixed Effects","Day-of-Month 5th Order Spline","Month Fixed Effects","Month x Year Effects","Week of Year FE"),covariate.labels=c(NA,"Mil. Pay Day/Day After","New Years Day/Eve", "Super Bowl", "Iditarod", "St Patricks Day", "Cinco de Mayo", "July 4th", "Labor Day Weekend", "Columbus Day Weekend", "Halloween and Weekend", "Thanksgiving", "Christmas", "Federal Holiday"), column.labels = names(FormulaList_Enforcement),no.space=T,keep.stat=c("n","adj.rsq"),notes="Newey-West Robust Errors.",type="text") 




rhs(FormulaList_Enforcement)=lapply(FormulaList_Enforcement,function(x) quote(WeekAftereWindow+MilitaryPayDay_ORdayAfter + 
                                                                            Holiday_NewYears + Holiday_SuperBowl + Holiday_Iditarod + 
                                                                            Holiday_StPat + Holiday_CincoDeMayo + Holiday_JulyForth + 
                                                                            Holiday_LaborOrWnknd + Holiday_ColumbusOrWnknd + Holiday_HalloweenOrWkndBF + 
                                                                            Holiday_Thanksgiving + Holiday_Christmas + Holiday_IsFederalHoliday + 
                                                                            poly(PRCP, 3) + poly(SNWD, 3) + poly(TMAX, 3) + DayOfWeek + 
                                                                            poly(dayofmonth, 5) + factor(MonthFixedEffect) * factor(YearFixedEffect)))

ModelList_enforcement=lapply(FormulaList_Enforcement,function(x) lm(x,subset(crimeCountsWide3)))
# Apply cluster robust method
DayAfterModelErrors_Enforce_VCOV=lapply(ModelList_enforcement, function(x) NeweyWest(x)) 
ClstrdErr_Enforce=lapply(DayAfterModelErrors_Enforce_VCOV, function(x) sqrt(diag(x)) )

##### *****Appendix C: Table C:2 ######
stargazer(ModelList_enforcement,omit=c("YearFixedEffect","poly","DayOfWeek","dayofmonth","MonthFixedEffect","factor.MonthFixedEffect...:factor.YearFixedEffect.","weekofyear"),omit.labels=c("Year Fixed Effects","Weather","Day of Week Fixed Effects","Day-of-Month 5th Order Spline","Month Fixed Effects","Month x Year Effects","Week of Year FE"),covariate.labels=c(NA,"Mil. Pay Day/Day After","New Years Day/Eve", "Super Bowl", "Iditarod", "St Patricks Day", "Cinco de Mayo", "July 4th", "Labor Day Weekend", "Columbus Day Weekend", "Halloween and Weekend", "Thanksgiving", "Christmas", "Federal Holiday"), column.labels = names(FormulaList_Enforcement),se=ClstrdErr_Enforce,no.space=T,keep.stat=c("n","adj.rsq"),notes="Newey-West Robust Errors.",type="text") 





########### **Results: Persistence ########### 


# Define time bins
crimeCountsWide3$num6weeksBefore=as.numeric(crimeCountsWide3$DaystoPFD<(43) & crimeCountsWide3$DaystoPFD>35)
crimeCountsWide3$num5weeksBefore=as.numeric(crimeCountsWide3$DaystoPFD<(36) & crimeCountsWide3$DaystoPFD>28)
crimeCountsWide3$num4weeksBefore=as.numeric(crimeCountsWide3$DaystoPFD<(29) & crimeCountsWide3$DaystoPFD>21)
crimeCountsWide3$num3weeksBefore=as.numeric(crimeCountsWide3$DaystoPFD<(22) & crimeCountsWide3$DaystoPFD>14)
crimeCountsWide3$num2weeksBefore=as.numeric(crimeCountsWide3$DaystoPFD<(15) & crimeCountsWide3$DaystoPFD>7)
crimeCountsWide3$num1weekBefore=as.numeric(crimeCountsWide3$DaystoPFD<(8) & crimeCountsWide3$DaystoPFD>0)

crimeCountsWide3$NumWeekAftereWindow=as.numeric(crimeCountsWide3$DaystoPFD<(0) & crimeCountsWide3$DaystoPFD>-8)
crimeCountsWide3$NumDay8to14=as.numeric(crimeCountsWide3$DaystoPFD<(-7) & crimeCountsWide3$DaystoPFD>-15)
crimeCountsWide3$NumDay15to21=as.numeric(crimeCountsWide3$DaystoPFD<(-14) & crimeCountsWide3$DaystoPFD>-22)
crimeCountsWide3$NumDay22to28=as.numeric(crimeCountsWide3$DaystoPFD<(-21) & crimeCountsWide3$DaystoPFD>-29)
crimeCountsWide3$NumDay29to36=as.numeric(crimeCountsWide3$DaystoPFD<(-28) & crimeCountsWide3$DaystoPFD>-36)



DesiredOutcomes=c("Violence","SubstanceFull","Property","Party","MedicalAssist")
EventVars=c("num4weeksBefore","num3weeksBefore","num2weeksBefore","num1weekBefore","NumWeekAftereWindow","NumDay8to14","NumDay15to21","NumDay22to28")


rhs(FormulaList)=lapply(FormulaList,function(x) quote(num4weeksBefore+num3weeksBefore+num2weeksBefore+num1weekBefore+NumWeekAftereWindow+NumDay8to14+NumDay15to21+NumDay22to28+MilitaryPayDay_ORdayAfter+Holiday_NewYears+Holiday_SuperBowl+Holiday_Iditarod+Holiday_StPat+Holiday_CincoDeMayo+Holiday_JulyForth+Holiday_LaborOrWnknd+Holiday_ColumbusOrWnknd+Holiday_HalloweenOrWkndBF+Holiday_Thanksgiving+Holiday_Christmas+Holiday_IsFederalHoliday+poly(PRCP,3)+poly(SNWD,3)+poly(TMAX,3)+DayOfWeek + poly(dayofmonth,5) +poly(as.numeric(dayofyear),5)*factor(YearFixedEffect) ))

FormulaList_Persistence_withMoPolyFEs=FormulaList[DesiredOutcomes]

ModelList=lapply(FormulaList_Persistence_withMoPolyFEs,function(x) lm(x,subset(crimeCountsWide3,Year>2009 & DaystoPFD!=0)))
ClstrdErrVCOV=lapply(ModelList, function(x) NeweyWest(x) )
ClstrdErr=lapply(ClstrdErrVCOV, function(x) sqrt(diag(x)) ) 

Persistance_Model_Results=ModelList
Persistance_Model_ClstrdErrVCOV=ClstrdErrVCOV
Persistance_Model_ClstrdErr=ClstrdErr


PersitanceModels=data.frame()
for(i in 1:length(Persistance_Model_Results)){
  PersitanceModels=rbind(PersitanceModels,data.frame(Outcome=names(Persistance_Model_Results)[i],coef(Persistance_Model_Results[[i]]),tidy(coefci(Persistance_Model_Results[[i]],vcov=Persistance_Model_ClstrdErrVCOV[[i]]))))
}

PersitanceModels=PersitanceModels[is.element(PersitanceModels$`.rownames`,EventVars),]
names(PersitanceModels)=c("Outcome","estimate","coef","x2.5","x97.5")
PersitanceModels$coef=factor(PersitanceModels$coef,levels = EventVars)

# Rename factor levels for plot
library(plyr)
PersitanceModels$coef=revalue(PersitanceModels$coef, c("num4weeksBefore"="-4","num3weeksBefore"="-3","num2weeksBefore"="-2","num1weekBefore"="-1","NumWeekAftereWindow"="PFD", "NumDay8to14"="1", "NumDay15to21"="2", "NumDay22to28"="3"))
PersitanceModels$Outcome=revalue(PersitanceModels$Outcome, c("SubstanceFull"="Substance", "MedicalAssist"="Medical Assist"))
detach("package:plyr", unload=TRUE)

PersitanceModels_Post2009_sample=PersitanceModels

##### *****Figure 1 ######
Figure1=ggplot(subset(PersitanceModels_Post2009_sample ),aes(x=coef,group=Outcome))+
  geom_hline(yintercept = 0,linetype=2)+geom_point(aes(y=estimate))+
  geom_line(aes(y=estimate))+ geom_errorbar(aes(ymin=x2.5, ymax= x97.5), width=.1)+
  ylab("Change in daily incidents and 95% Conf.") + 
  coord_cartesian(ylim=c(-7,15)) + 
  facet_wrap(~Outcome)+
  theme_bw()+
  xlab("")

ggsave("./Figures/Figure1.pdf",plot=Figure1, height=4, width=6, units='in')


# Full Sample persistence

ModelList=lapply(FormulaList_Persistence_withMoPolyFEs,function(x) lm(x,subset(crimeCountsWide3, DaystoPFD!=0)))
ClstrdErrVCOV=lapply(ModelList, function(x) NeweyWest(x) )
ClstrdErr=lapply(ClstrdErrVCOV, function(x) sqrt(diag(x)) ) 

Persistance_Model_Results_FullSampe=ModelList
Persistance_Model_ClstrdErrVCOV_FullSampe=ClstrdErrVCOV
Persistance_Model_ClstrdErr_FullSampe=ClstrdErr

# Full sample table:
stargazer(Persistance_Model_Results_FullSampe,omit=c("YearFixedEffect","poly","Holiday_","DayOfWeek","dayofmonth","MonthFixedEffect","factor.MonthFixedEffect...:factor.YearFixedEffect."),omit.labels=c("Year Fixed Effects","Weather","Holiday Effects","Day of Week Fixed Effects","Day-of-Month 5th Order Spline","Month Fixed Effects","Month x Year Effects"), column.labels = names(Persistance_Model_Results_FullSampe),se=Persistance_Model_ClstrdErr_FullSampe,no.space=T,keep.stat=c("n","adj.rsq"),notes="Newey-West Robust Errors.") 


####### Plot of persistence in full sample  

# PersitanceModels=data.frame()
# for(i in 1:length(Persistance_Model_Results)){
#   PersitanceModels=rbind(PersitanceModels,data.frame(Outcome=names(Persistance_Model_Results_FullSampe)[i],coef(Persistance_Model_Results_FullSampe[[i]]),tidy(coefci(Persistance_Model_Results_FullSampe[[i]],vcov=Persistance_Model_ClstrdErrVCOV_FullSampe[[i]]))))
# }
# 
# PersitanceModels=PersitanceModels[is.element(PersitanceModels$`.rownames`,EventVars),]
# names(PersitanceModels)=c("Outcome","estimate","coef","x2.5","x97.5")
# PersitanceModels$coef=factor(PersitanceModels$coef,levels = EventVars)
# 
# library(plyr)
# PersitanceModels$coef=revalue(PersitanceModels$coef, c("num4weeksBefore"="-4","num3weeksBefore"="-3","num2weeksBefore"="-2","num1weekBefore"="-1","NumWeekAftereWindow"="PFD", "NumDay8to14"="1", "NumDay15to21"="2", "NumDay22to28"="3"))
# PersitanceModels$Outcome=revalue(PersitanceModels$Outcome, c("SubstanceFull"="Substance", "MedicalAssist"="Medical Assist"))
# detach("package:plyr", unload=TRUE)
# 
# PersitanceModels_Full_sample=PersitanceModels
# PersitanceModels_Full_sample$Payments="All Years"
# 
# # Plot of persistence in full sample
# ggplot(subset(PersitanceModels_Full_sample ),aes(x=coef,group=Outcome))+geom_hline(yintercept = 0,linetype=2)+geom_point(aes(y=estimate))+geom_line(aes(y=estimate))+ geom_errorbar(aes(ymin=x2.5, ymax= x97.5), width=.1)+ylab("Change in daily call rate and 95% Conf.") + coord_cartesian(ylim=c(-7,15)) + facet_wrap(~Outcome)+theme_bw()+xlab("")
# 
#
# # # # # # # # # # # # # # ## # # # # ## # # # # ## # # # # # # # 


# Full and subsetted sample tables:
##### *****Appendix: Table 6 ######

# Top part of table:

stargazer(list(Persistance_Model_Results,Persistance_Model_Results_FullSampe),omit=c("YearFixedEffect","poly","Holiday_","DayOfWeek","dayofmonth","MonthFixedEffect","factor.MonthFixedEffect...:factor.YearFixedEffect."), omit.labels=c("Year Fixed Effects","Weather","Holiday Effects","Day of Week Fixed Effects","Day-of-Month 5th Order Spline","Month Fixed Effects","Month x Year Effects"), column.labels = names(c(Persistance_Model_Results,Persistance_Model_Results_FullSampe)),se=c(Persistance_Model_ClstrdErr,Persistance_Model_ClstrdErr_FullSampe),no.space=T,keep.stat=c("n","adj.rsq"),notes="Newey-West Robust Errors.") 

# Bottom part of table:

# Set hypotheses tests for the cumulative effects
hypoth_2="7*NumWeekAftereWindow+7*NumDay8to14=0"
hypoth_3="7*NumWeekAftereWindow+7*NumDay8to14+7*NumDay15to21=0"
hypoth_4="7*NumWeekAftereWindow+7*NumDay8to14+7*NumDay15to21+7*NumDay22to28=0"

# Print estimates and 95% ci's for the table at 2, 3, and 4 weeks
cumweeks_2=mapply(function(X,Y) as.data.frame(tidy(confint((glht(X, linfct = hypoth_2, vcov =Y ))))[,c("estimate","conf.low","conf.high")]), X=c(Persistance_Model_Results,Persistance_Model_Results_FullSampe),Y=c(Persistance_Model_ClstrdErrVCOV,Persistance_Model_ClstrdErrVCOV_FullSampe))

cumweeks_2_est=sprintf("%.2f",unlist(lapply(as.data.frame(cumweeks_2), function(x) x$estimate[1])))
cumweeks_2_ci=paste("[",sprintf("%.1f",unlist(lapply(as.data.frame(cumweeks_2), function(x) x$conf.low[1]))),",",sprintf("%.1f",unlist(lapply(as.data.frame(cumweeks_2), function(x) x$conf.high[1]))),"]",sep="")


cumweeks_3=mapply(function(X,Y) as.data.frame(tidy(confint((glht(X, linfct = hypoth_3, vcov =Y ))))[,c("estimate","conf.low","conf.high")]), X=c(Persistance_Model_Results,Persistance_Model_Results_FullSampe),Y=c(Persistance_Model_ClstrdErrVCOV,Persistance_Model_ClstrdErrVCOV_FullSampe))

cumweeks_3_est=sprintf("%.2f",unlist(lapply(as.data.frame(cumweeks_3), function(x) x$estimate[1])))
cumweeks_3_ci=paste("[",sprintf("%.1f",unlist(lapply(as.data.frame(cumweeks_3), function(x) x$conf.low[1]))),",",sprintf("%.1f",unlist(lapply(as.data.frame(cumweeks_3), function(x) x$conf.high[1]))),"]",sep="")


cumweeks_4=mapply(function(X,Y) as.data.frame(tidy(confint((glht(X, linfct = hypoth_4, vcov =Y ))))[,c("estimate","conf.low","conf.high")]), X=c(Persistance_Model_Results,Persistance_Model_Results_FullSampe),Y=c(Persistance_Model_ClstrdErrVCOV,Persistance_Model_ClstrdErrVCOV_FullSampe))

cumweeks_4_est=sprintf("%.2f",unlist(lapply(as.data.frame(cumweeks_4), function(x) x$estimate[1])))
cumweeks_4_ci=paste("[",sprintf("%.1f",unlist(lapply(as.data.frame(cumweeks_4), function(x) x$conf.low[1]))),",",sprintf("%.1f",unlist(lapply(as.data.frame(cumweeks_4), function(x) x$conf.high[1]))),"]",sep="")


# Latex Cumulative effects table:

print(xtable(rbind(
  data.frame(variable="Total Cumulative Effect at Week 2:",t(cumweeks_2_est)) ,
  data.frame(variable="95% Conf.", t(cumweeks_2_ci)),

  data.frame(variable="Total Cumulative Effect at Week 3:",t(cumweeks_3_est)) ,
  data.frame(variable="95% Conf.", t(cumweeks_3_ci)),
  
  data.frame(variable="Total Cumulative Effect at Week 4:",t(cumweeks_4_est)) ,
  data.frame(variable="95% Conf.", t(cumweeks_4_ci))
  
  )), include.rownames=FALSE)

write.table(rbind(
  data.frame(variable="Total Cumulative Effect at Week 2:",t(cumweeks_2_est)) ,
  data.frame(variable="95% Conf.", t(cumweeks_2_ci)),
  
  data.frame(variable="Total Cumulative Effect at Week 3:",t(cumweeks_3_est)) ,
  data.frame(variable="95% Conf.", t(cumweeks_3_ci)),
  
  data.frame(variable="Total Cumulative Effect at Week 4:",t(cumweeks_4_est)) ,
  data.frame(variable="95% Conf.", t(cumweeks_4_ci))
  
), "clipboard", sep="\t", row.names=FALSE)




# Joint significance test

Hnull <- c("NumWeekAftereWindow=0","NumDay8to14=0","NumDay15to21=0","NumDay22to28=0")

JointFtest=mapply(function(X,Y) data.frame(Fstat=linearHypothesis(X,Hnull,vcov=Y)$F[2],Pval=linearHypothesis(X,Hnull,vcov=Y)$`Pr(>F)`[2]),X=c(Persistance_Model_Results,Persistance_Model_Results_FullSampe),Y=c(Persistance_Model_ClstrdErrVCOV,Persistance_Model_ClstrdErrVCOV_FullSampe))

JointFtest=as.data.frame(t(as.data.frame(JointFtest)))

JointFtest$Stars=""
JointFtest$Stars[JointFtest$Pval<.1]="*"
JointFtest$Stars[JointFtest$Pval<.05]="**"
JointFtest$Stars[JointFtest$Pval<.01]="***"
JointFtest$Fstat=paste(sprintf("%.3f",JointFtest$Fstat),JointFtest$Stars,sep="")
JointFtest$'Unadjusted P-value'=paste("[",sprintf("%.3f",JointFtest$Pval),"]",sep="")
JointFtest$'Adjusted P-value'=paste("[",sprintf("%.3f",p.adjust(JointFtest$Pval)),"]",sep="")


xtable(t(JointFtest[,c("Fstat", 'Unadjusted P-value', 'Adjusted P-value')]))



#Economic Signficance. Calculating the % change over the annual average for the 4 weeks after distribution
TotalCrimebyCat2016=data.frame(round(t(unlist(lapply(FormulaList[DesiredOutcomes], function(x) sum(rowSums(subset(crimeCountsWide3,Year==2016)[,lhs.vars(x),drop=FALSE])) )) ),2))

Pc_Point_cm=unlist(c(unlist(lapply(as.data.frame(cumweeks_4), function(x) x$estimate[1]))[1:5]*100/TotalCrimebyCat2016,unlist(lapply(as.data.frame(cumweeks_4), function(x) x$estimate[1]))[6:10]*100/TotalCrimebyCat2016))
Pc_Low_cm=unlist(c(unlist(lapply(as.data.frame(cumweeks_4), function(x) x$conf.low[1]))[1:5]*100/TotalCrimebyCat2016,unlist(lapply(as.data.frame(cumweeks_4), function(x) x$conf.low[1]))[6:10]*100/TotalCrimebyCat2016))
Pc_High_cm=unlist(c(unlist(lapply(as.data.frame(cumweeks_4), function(x) x$conf.high[1]))[1:5]*100/TotalCrimebyCat2016,unlist(lapply(as.data.frame(cumweeks_4), function(x) x$conf.high[1]))[6:10]*100/TotalCrimebyCat2016))

rbind(Pc_Low_cm,Pc_Point_cm,Pc_High_cm)

paste("[",sprintf("%.2f",Pc_Low_cm),"%","]",sep="")
paste(sprintf("%.2f",Pc_Point_cm),"%",sep="")

PcChangTab=data.frame(rbind(TotalCrimebyCat2016,sprintf("%.2f",cumweeks_4["estimate",1:5]),
paste(sprintf("%.2f",Pc_Point_cm),"%",sep=""),
paste("[",sprintf("%.2f",Pc_Low_cm),"%",",",sprintf("%.2f",Pc_High_cm),"%","]",sep="")))

#names(PcChangTab)=names(Pc_Point_cm)

PcChangTab[,1:5]

row.names(PcChangTab)=c("Annual Incidents, 2016","Cummulative Change After PFD","% Change After PFD","95% Confidence")
xtable(PcChangTab)




########### **Results: Daily persistance effect  ########### 


crimeCountsWide3$PFDDaily="AAllOtherDays"
crimeCountsWide3$PFDDaily[which(abs(crimeCountsWide3$DaystoPFD)<=30)]=as.character(crimeCountsWide3$DaystoPFD[which(abs(crimeCountsWide3$DaystoPFD)<=30)])


crimeCountsWide3$PFDDaily[which(crimeCountsWide3$DaystoPFD>=-28 & crimeCountsWide3$DaystoPFD<0)]=as.character(crimeCountsWide3$DaystoPFD[which(crimeCountsWide3$DaystoPFD>=-28 & crimeCountsWide3$DaystoPFD<0)])
crimeCountsWide3$PFDDaily[which(crimeCountsWide3$DaystoPFD<=28 & crimeCountsWide3$DaystoPFD>=0)]=as.character(crimeCountsWide3$DaystoPFD[which(crimeCountsWide3$DaystoPFD<=28 & crimeCountsWide3$DaystoPFD>=0)])
LVLS=levels(factor(crimeCountsWide3$PFDDaily))
crimeCountsWide3$PFDDaily=factor(crimeCountsWide3$PFDDaily,levels=sort(LVLS,decreasing = T))

FormulaList_dailyevent=FormulaList[DesiredOutcomes]

rhs(FormulaList_dailyevent)=lapply(FormulaList_dailyevent,function(x) quote(PFDDaily+MilitaryPayDay_ORdayAfter+MilitaryPayDay_ORdayAfter+Holiday_NewYears+Holiday_SuperBowl+Holiday_Iditarod+Holiday_StPat+Holiday_CincoDeMayo+Holiday_JulyForth+Holiday_LaborOrWnknd+Holiday_ColumbusOrWnknd+Holiday_HalloweenOrWkndBF+Holiday_Thanksgiving+Holiday_Christmas+Holiday_IsFederalHoliday+poly(PRCP,3)+poly(SNWD,3)+poly(TMAX,3)+DayOfWeek + poly(dayofmonth,5) +poly(as.numeric(dayofyear),5)*factor(YearFixedEffect)))



ModelList=lapply(FormulaList_dailyevent,function(x) lm(x,subset(crimeCountsWide3,Year>=2009 & DaystoPFD!=0)))
ClstrdErrVCOV=lapply(ModelList, function(x) NeweyWest(x) )

ModelList_dailyevent=lapply(ModelList,function(x) tidy(x,conf.int = T,level=.9))
CoefList_dailyevent=lapply(ModelList_dailyevent,function(x) subset(x,grepl("PFDDaily",x$term)))

CoefList_dailyevent_a=data.frame()
for(i in 1:length(CoefList_dailyevent)){
  CoefList_dailyevent[[i]]$term=as.numeric(gsub("PFDDaily","",CoefList_dailyevent[[i]]$term))
  CoefList_dailyevent_a=rbind(CoefList_dailyevent_a,data.frame(Outcome=names(CoefList_dailyevent)[[i]],CoefList_dailyevent[[i]]))
}

CoefList_dailyevent_a=CoefList_dailyevent_a[order(CoefList_dailyevent_a$Outcome,CoefList_dailyevent_a$term),]
write.csv(CoefList_dailyevent_a,"C:/Users/brett/Google Drive/Daily Crime Data (2012 to 2016)/DailyChange.csv")

WeekendColors=unique(crimeCountsWide3$PFDDaily[which(abs(crimeCountsWide3$DaystoPFD)<=30 & crimeCountsWide3$DayOfWeek %in% c("Saturday","Sunday") & crimeCountsWide3$Year>2009)])




Daily_PersitanceModels=data.frame()
for(i in 1:length(Persistance_Model_Results)){
  Daily_PersitanceModels=rbind(Daily_PersitanceModels,data.frame(Outcome=names(ModelList)[i],coef(ModelList[[i]]),tidy(coefci(ModelList[[i]],vcov=ClstrdErrVCOV[[i]]))))
}

Daily_PersitanceModels=subset(Daily_PersitanceModels,grepl("PFDDaily", .rownames))
Daily_PersitanceModels$.rownames=as.numeric(gsub("PFDDaily","",Daily_PersitanceModels$.rownames))
Daily_PersitanceModels=Daily_PersitanceModels[order(Daily_PersitanceModels$Outcome,Daily_PersitanceModels$.rownames),]

##### *****Appendix: Figure A:2 ######
FigureA2=ggplot(Daily_PersitanceModels,aes(y=coef.ModelList..i...,x=-.rownames))+
  geom_point(aes(group=1,col=.rownames%in%WeekendColors),size = 1)+geom_hline(yintercept=0)+
  scale_fill_brewer(name="")+
  geom_errorbar(aes(ymin=X2.5..,ymax=X97.5..,col=.rownames%in%WeekendColors),size = 0.25)+
  geom_smooth(aes(group=.rownames>0,fill="Local smoother"),col="darkgrey",size = 1,se=F)+
  facet_wrap(~Outcome)+
  theme_bw()+ 
  scale_colour_Publication(name = "Day of Week",labels=c("Weekday","Weekend"))+
  ylab("Change in daily incidents and 95% confin")+xlab("Days Before/After PFD Payment")+
  coord_cartesian(ylim=c(-15,20))+
  theme(legend.position = c(.8,.2))

FigureA2
ggsave("./Figures/FigureA2.pdf",plot=FigureA2, height=4, width=6.5, units='in')



########### **Results: Marginal Effect  ########### 


DesiredOutcomes=c("Violence","SubstanceFull","Property","Party","MedicalAssist")

rhs(FormulaList)=lapply(FormulaList,function(x) quote(MilitaryPayDay_ORdayAfter+Holiday_NewYears+Holiday_SuperBowl+Holiday_Iditarod+Holiday_StPat+Holiday_CincoDeMayo+Holiday_JulyForth+Holiday_LaborOrWnknd+Holiday_ColumbusOrWnknd+Holiday_HalloweenOrWkndBF+Holiday_Thanksgiving+Holiday_Christmas+Holiday_IsFederalHoliday+poly(PRCP,3)+poly(SNWD,3)+poly(TMAX,3)+DayOfWeek + poly(dayofmonth,5) + factor(MonthFixedEffect)*factor(YearFixedEffect)))

#NumPFDdayAfter+NumPFDdayAfter:Amount_First_Date_100m+NumPFDdayAfter:I(Amount_First_Date_100m^2)

PersistanceWindow=c("NumPFDdayAfter","NumPFDweekend","NumWeekAftereWindow")
PersistanceWindow=c("NumWeekAftereWindow")


FormulaList2=FormulaList[DesiredOutcomes]
PersistanceWindow_Models=list()
PersistanceWindow_ClstrdErrVCOV=list()
PersistanceWindow_ClstrdErr=list()
Ouctomes_byWindow_atmeans=data.frame()
Ouctomes_byWindow_atobs=data.frame()

for(i in 1:length(PersistanceWindow)){

    treatment=paste(paste(PersistanceWindow[i],c("",":Amount_First_Date_100m",":I(Amount_First_Date_100m^2)"),sep=""),collapse = " + ")
  rhs(FormulaList2)=lapply(FormulaList[DesiredOutcomes],function(x) rhs(as.formula(paste("~",treatment,"+",paste(as.character(rhs(FormulaList[[1]]))[-1],collapse="+") ))) )
  
  ModelList=lapply(FormulaList2,function(x) lm(x,subset(crimeCountsWide3, DaystoPFD!=0 )))
  ClstrdErrVCOV=lapply(ModelList, function(x) NeweyWest(x) )
  MrginalFX_Persist=data.frame()
  
  pfdpayouts=unique(subset(crimeCountsWide3)$Amount_First_Date_100m)
  hypoth=paste("1*",PersistanceWindow[i],":Amount_First_Date_100m + ",2*pfdpayouts,"*",PersistanceWindow[i],":I(Amount_First_Date_100m^2)=0",sep="")

  for(j in 1:length(names(FormulaList2))){
    MrginalFX_i=data.frame(do.call(rbind, lapply(lapply(hypoth,function(x) tidy(confint(glht(ModelList[[j]], linfct = x, vcov = ClstrdErrVCOV[[j]] ), level = 0.90))), data.frame, stringsAsFactors=FALSE)),PFD_Payouts=pfdpayouts)
    MrginalFX_Persist=rbind(MrginalFX_Persist,data.frame(Outcome=names(FormulaList)[j],Window=PersistanceWindow[i],MrginalFX_i))
  }
  Ouctomes_byWindow_atobs=rbind(Ouctomes_byWindow_atobs,MrginalFX_Persist)
  
  MrginalFX_Persist_atmean=data.frame() 
  hypoth=paste("1*",PersistanceWindow[i],":Amount_First_Date_100m + ",2*mean(unique(subset(crimeCountsWide3)$Amount_First_Date_100m)),"*",PersistanceWindow[i],":I(Amount_First_Date_100m^2)=0",sep="")
  for(j in 1:length(names(FormulaList2))){
    MrginalFX_i= tidy(summary(glht(ModelList[[j]], linfct = hypoth, vcov = ClstrdErrVCOV[[j]] )))
    MrginalFX_Persist_atmean=rbind(MrginalFX_Persist_atmean,data.frame(Outcome=names(FormulaList2)[j],Window=PersistanceWindow[i],MrginalFX_i))
  }
  
  Ouctomes_coefs=MrginalFX_Persist_atmean
  
  Ouctomes_coefs$Stars=""
  Ouctomes_coefs$Stars[Ouctomes_coefs$'p.value'<.1]="*"
  Ouctomes_coefs$Stars[Ouctomes_coefs$'p.value'<.05]="**"
  Ouctomes_coefs$Stars[Ouctomes_coefs$'p.value'<.01]="***"
  Ouctomes_coefs$EstimateSig=paste(round(Ouctomes_coefs$estimate,3),Ouctomes_coefs$Stars,sep="")
  Ouctomes_coefs$PrintError=paste("(",round(Ouctomes_coefs$"std.error",3),")",sep="")
  Ouctomes_coefs$UnadjustedP=paste("[",format(round(Ouctomes_coefs$"p.value",3),digits=2),"]",sep="")
  Ouctomes_coefs$adjustP=paste("[",format(round(p.adjust(Ouctomes_coefs$"p.value"),3),digits=2),"]",sep="")
  Ouctomes_coefsL=melt(Ouctomes_coefs[,c("Outcome","EstimateSig","PrintError","UnadjustedP","adjustP")],id.vars="Outcome")
  Ouctomes_coefsL$Outcome=factor(Ouctomes_coefsL$Outcome,levels=names(ModelList))
  Ouctomes_coefsW=dcast(Ouctomes_coefsL,variable~Outcome)
  Ouctomes_coefsW$variable=as.character(Ouctomes_coefsW$variable)
  Ouctomes_coefsW[1,1]=PersistanceWindow[i]
  Ouctomes_byWindow_atmeans=rbind(Ouctomes_byWindow_atmeans,Ouctomes_coefsW)
  
  #PersistanceWindow_Models[[i]]=ModelList
  #PersistanceWindow_ClstrdErrVCOV[[i]]=PersistanceWindow_ClstrdErrVCOV
  #PersistanceWindow_ClstrdErr[[i]]=PersistanceWindow_ClstrdErr
  
}

MarginalMainPersistenceEffects_atmeans_FullSample=Ouctomes_byWindow_atmeans
MarginalMainPersistenceEffects_FullSample=Ouctomes_byWindow_atobs


###########  *****Table 3  ########### 
print(xtable(MarginalMainPersistenceEffects_atmeans_FullSample), include.rownames=FALSE)
print(xtable(rbind(
  data.frame(variable="Mean Daily Call Count",round(t(unlist(lapply(FormulaList[DesiredOutcomes], function(x) mean(rowSums(subset(crimeCountsWide3)[,lhs.vars(x),drop=FALSE])) )) ),2)),
  data.frame(variable="St. dev Call Count",round(t(unlist(lapply(FormulaList[DesiredOutcomes], function(x) sd(rowSums(subset(crimeCountsWide3)[,lhs.vars(x),drop=FALSE])) ))),2)))))






########### **Results: Compare to SNAP ########### 
FormulaList_log_SNAP=FormulaList
crimeCountsWide4=crimeCountsWide3

crimeCountsWide4$Month3=substr(format(as.Date(crimeCountsWide4$Received_Date),format="%B"),1,3)
crimeCountsWide4=merge(crimeCountsWide4,SnapPayments,by.x=c("Year","Month3"),by.y=c("Year","Month"))

crimeCountsWide4$SNAP_numeric=as.numeric(crimeCountsWide4$dayofmonth<8)


# Joitly Estimate PFD and SNAP marginal effect 
DesiredOutcomes=c("Violence","SubstanceFull","Property","Party","MedicalAssist")
rhs(FormulaList)=lapply(FormulaList[DesiredOutcomes],function(x) quote(NumWeekAftereWindow+NumWeekAftereWindow:I(PFD_Amount_FromAvg*100)+SNAP_numeric+SNAP_numeric:I((CostDec2016-11000366)/1000000)+MilitaryPayDay_ORdayAfter+Holiday_NewYears+Holiday_SuperBowl+Holiday_Iditarod+Holiday_StPat+Holiday_CincoDeMayo+Holiday_JulyForth+Holiday_LaborOrWnknd+Holiday_ColumbusOrWnknd+Holiday_HalloweenOrWkndBF+Holiday_Thanksgiving+Holiday_Christmas+Holiday_IsFederalHoliday+poly(PRCP,3)+poly(SNWD,3)+poly(TMAX,3)+DayOfWeek +factor(MonthFixedEffect)*factor(YearFixedEffect)))


ModelList=lapply(FormulaList[DesiredOutcomes],function(x) lm(x,subset(crimeCountsWide4, DaystoPFD!=0)))
ClstrdErrVCOV=lapply(ModelList, function(x) NeweyWest(x) )
ClstrdErr=lapply(ClstrdErrVCOV, function(x) sqrt(diag(x)) )

###########  *****Table 3  ########### 
# top of table
stargazer(ModelList,omit=c("YearFixedEffect","poly","Holiday","DayOfWeek","dayofmonth","MonthFixedEffect","factor.MonthFixedEffect...:factor.YearFixedEffect."),omit.labels=c("Year Fixed Effects","Weather","Holiday","Day of Week Fixed Effects","Day-of-Month 5th Order Spline","Month Fixed Effects","Month x Year Effects"),se=ClstrdErr, column.labels = DesiredOutcomes,no.space=T,keep.stat=c("n","adj.rsq"),notes="Newey-West Robust Errors.",type="text") 


# top of table,  Elasticities

# Mean values
DailyCrimeMeans=1/round(unlist(lapply(FormulaList[DesiredOutcomes], function(x) mean(rowSums(subset(crimeCountsWide3)[,lhs.vars(x),drop=FALSE])) )),2)

# % change in dollars at mean
1/697 # average pfd distribution (rounded)
1/10 # average SNAP distribution (rounded)
Scalers=data.frame(PFD_Scale=DailyCrimeMeans/(1/697),SNAP_SCALE=DailyCrimeMeans/(1/11))

hypoth="meanIncident_SNAP*SNAP_numeric:I((CostDec2016 - 11000366)/1e+06) - meanIncident_PFD*NumWeekAftereWindow:I(PFD_Amount_FromAvg * 100)=0"
hypoth=apply(data.frame(Scalers,hypoth),1,function(x) gsub("meanIncident_PFD",x["PFD_Scale"], gsub("meanIncident_SNAP",x["SNAP_SCALE"],x["hypoth"] )) )

SNAPvPFDtest=mapply(function(X,Y,Z) tidy(summary(glht(X, linfct = Z, vcov = Y)))[,-1],X=ModelList,Y=ClstrdErrVCOV,Z=hypoth)

Elasticities=mapply(function(X,Y,Z) data.frame(E_SNAP=coef(X)["SNAP_numeric:I((CostDec2016 - 11000366)/1e+06)"]*Z,E_PFD=coef(X)["NumWeekAftereWindow:I(PFD_Amount_FromAvg * 100)"]*Y ),X=ModelList,Y=Scalers$PFD_Scale,Z=Scalers$SNAP_SCALE)

xtable(rbind(Elasticities,SNAPvPFDtest),digits=3)
 # Bonfironi corrections:
xtable(t(data.frame(AdjP=paste("[",format(round(p.adjust(t(SNAPvPFDtest)[,"p.value"]),3),digits=2),"]",sep=""))))


##### Placebo #####

#DesiredOutcomes=c("Violence","SubstanceFull","Property","Party","MedicalAssist")

rhs(FormulaList)=lapply(FormulaList,function(x) quote(I(WeeksToPFD==i)++MilitaryPayDay_ORdayAfter+Holiday_NewYears+Holiday_SuperBowl+Holiday_Iditarod+Holiday_StPat+Holiday_CincoDeMayo+Holiday_JulyForth+Holiday_LaborOrWnknd+Holiday_ColumbusOrWnknd+Holiday_HalloweenOrWkndBF+Holiday_Thanksgiving+Holiday_Christmas+Holiday_IsFederalHoliday+poly(PRCP,3)+poly(SNWD,3)+poly(TMAX,3)+DayOfWeek + poly(dayofmonth,5) +factor(MonthFixedEffect)*factor(YearFixedEffect)))

crimeCountsWide3$WeeksToPFD=as.numeric(ceiling((crimeCountsWide3$DaystoPFD+1)/7))
#table(crimeCountsWide3$WeeksToPFD)

PlaceboResults=data.frame()
for(i in min(crimeCountsWide3$WeeksToPFD):max(crimeCountsWide3$WeeksToPFD)){
  for(j in 1:length(FormulaList[DesiredOutcomes])){
    rsdtmp=tidy(lm(FormulaList[DesiredOutcomes][[j]],subset(crimeCountsWide3)), conf.int=T,conf.level = 0.95, func = lmtest::coefci, vcov = NeweyWest)
    rsdtmp$WeeksToPFD=i
    rsdtmp$Outcome=names(FormulaList[DesiredOutcomes])[[j]]
    PlaceboResults=rbind(PlaceboResults,rsdtmp)
  }
}
PlaceboResults=subset(PlaceboResults,Outcome%in%DesiredOutcomes)

PlaceboResults$Outcome=factor(PlaceboResults$Outcome,levels=c("Violence","SubstanceFull","Property","Party","MedicalAssist"))

gg = ggplot(subset(PlaceboResults,term=="I(WeeksToPFD == i)TRUE"),aes(y=estimate,x=-WeeksToPFD))+ geom_vline(xintercept=0)+geom_line()+geom_point(aes(alpha=p.value<.05))
gg = gg + geom_line(aes(y=`conf.high`,x=-WeeksToPFD),lty=2)
gg = gg + geom_line(aes(y=`conf.low`,x=-WeeksToPFD),lty=2)
gg = gg + geom_hline(yintercept=0)
gg = gg + facet_wrap(~Outcome)
gg = gg + coord_cartesian(ylim=c(-7.5,7.5))
gg = gg + scale_x_continuous(breaks=scales::pretty_breaks(n=5)) 
gg = gg + xlab("Weeks from PFD Distribution")
gg = gg + ylab("Estimate and 95% Conf.")
gg = gg + theme(legend.position = "none")
gg = gg + theme_bw()
gg
#Daylight savings?

PlaceboResultsTRUE=PlaceboResults[which(PlaceboResults$WeeksToPFD==0),c("term","Outcome","estimate")]
names(PlaceboResultsTRUE)[3]="TrueTreatmentEff"
PlaceboResultsTRUE=merge(PlaceboResults,PlaceboResultsTRUE,by=c("term","Outcome"))


subset(PlaceboResults,term=="I(WeeksToPFD == i)TRUE" & Outcome=="SubstancePart")[order(abs(subset(PlaceboResults,term=="I(WeeksToPFD == i)TRUE" & Outcome=="SubstancePart")$estimate)),]


PlaceboResultsTRUE=mutate(group_by(PlaceboResultsTRUE,term,Outcome),ProbabilityValue=mean(abs(estimate) >= abs(TrueTreatmentEff)))

###########  *****Figure 2  ########### 
gg = ggplot(subset(PlaceboResultsTRUE,term=="I(WeeksToPFD == i)TRUE"&Outcome!="TrafficParking"))+geom_histogram(aes(x=abs(estimate)))+ facet_wrap(~Outcome)
gg = gg + geom_vline(data=subset(PlaceboResultsTRUE,term=="I(WeeksToPFD == i)TRUE"&Outcome!="TrafficParking" & WeeksToPFD==0),aes(xintercept = abs(TrueTreatmentEff)))
gg = gg + xlab("Placebo and True Treatment Effect Magnitude")
gg = gg + theme_bw()
gg = gg + geom_vline(data=subset(PlaceboResultsTRUE,term=="I(WeeksToPFD == i)TRUE"&Outcome!="TrafficParking" & WeeksToPFD==0),aes(xintercept = abs(TrueTreatmentEff)))
gg = gg + geom_text(data=subset(PlaceboResultsTRUE,term=="I(WeeksToPFD == i)TRUE"&Outcome!="TrafficParking" & WeeksToPFD==0),aes(x = abs(TrueTreatmentEff)+1.1,y=15,label = paste("p=",round(ProbabilityValue,3))))
gg = gg +coord_cartesian(xlim=c(0,6))
Figure2=gg


ggsave("./Figures/Figure2.pdf",plot=Figure2, height=4, width=6.5, units='in')



########  Plot Call Data  ######### 



pfdDatesAmounts_checks2=data.frame()
for(i in 0:6){
  pfdDatesAmounts_checks2=rbind(pfdDatesAmounts_checks2,data.frame(pfdDatesAmounts_checks[,c(-1,-5)],Date=pfdDatesAmounts_checks$PFDSend_Date+1))
}

pfdDatesAmounts_checks2$Persons_zscr=(pfdDatesAmounts_checks2$Persons-mean(pfdDatesAmounts_checks2$Persons))/sd(pfdDatesAmounts_checks2$Persons)
pfdDatesAmounts_checks2$PFD_Amount_zscr=(pfdDatesAmounts_checks2$PFD_Amount-mean(pfdDatesAmounts_checks2$PFD_Amount))/sd(pfdDatesAmounts_checks2$PFD_Amount)
pfdDatesAmounts_checks2$Total_Amount_zscr=(pfdDatesAmounts_checks2$Total_Amount-mean(pfdDatesAmounts_checks2$Total_Amount))/sd(pfdDatesAmounts_checks2$Total_Amount)


crimeCountsWide5=merge(crimeCountsWide3,pfdDatesAmounts_checks2,by.x="Received_Date",by.y="Date",all.x=T)
crimeCountsWide5$PFD_Amount.y[is.na(crimeCountsWide5$PFD_Amount.y)]=0
crimeCountsWide5$Total_Amount[is.na(crimeCountsWide5$Total_Amount)]=0
crimeCountsWide5$Persons[is.na(crimeCountsWide5$Persons)]=0

crimeCountsWide5$PFD_Amount_zscr[is.na(crimeCountsWide5$PFD_Amount_zscr)]=0
crimeCountsWide5$Total_Amount_zscr[is.na(crimeCountsWide5$Total_Amount_zscr)]=0
crimeCountsWide5$Persons_zscr[is.na(crimeCountsWide5$Persons_zscr)]=0



AggCountsByOutnDay=lapply(FormulaList[DesiredOutcomes], function(x) data.frame(crimeCountsWide5[,c("Received_Date","Year","DayofYear","dayofyear","dayofmonth","DayOfWeekShort","Month","PFD_Amount.y")],Calls=rowSums(crimeCountsWide5[,lhs.vars(x),drop=FALSE])))
namesTitle=list("Violence","Substance","Property","Party","Medical Assist")


for(i in 1:length(namesTitle)){
  AggCountsByOutnDay[[i]]=data.frame(AggCountsByOutnDay[[i]],Outcome=namesTitle[[i]])
}


# Day of week plot
DOWplot=lapply(AggCountsByOutnDay,function(X) ggplot(X,aes(x=(DayOfWeekShort),y=Calls, group=1))+
                 stat_summary(fun.y = "median", geom = "line")+
                 xlab(NULL)+ylab("Calls/day")+
                 theme_bw()+ theme(plot.margin = margin(t=1,r=1,b=1,l=1))+ 
                 annotate(geom="text" ,x="Wed",y=Inf,label=unique(X$Outcome),vjust=1.2 )+ 
                 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) )

# month of year plot

MOYplot=lapply(AggCountsByOutnDay,function(X) ggplot(X,aes(x=Month,y=Calls,group=1))+stat_summary(fun.y = "median", geom = "line")+xlab(NULL)+ylab("Calls/day")+theme_bw()+ theme(plot.margin = margin(t=1,r=1,b=1,l=1))+ annotate(geom="text" ,x=07,y=Inf,label=unique(X$Outcome),vjust=1.2 )+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()))



g1 <- arrangeGrob(grobs=DOWplot, ncol=1, padding = unit(0, "line"))
g2 <- arrangeGrob(grobs=MOYplot, ncol=1, padding = unit(0, "line"))

###########  *****Appendix: Figure A:1  ########### 
grid.arrange(g1, g2, ncol=2)

ggsave("./Figures/FigureA1.pdf",plot=grid.arrange(g1, g2, ncol=2), height=6, width=6, units='in')



######## Call Codes and Categories ######### 

###########  *****Appendix: Table A:1  ########### 
CallAssignements=unlist(lapply(FormulaList, function(x) lhs.vars(x)))
CallAssignements2=data.frame(OutcomeCatagory=names(CallAssignements),CallCode=CallAssignements,`Average Daily Calls`=colMeans(crimeCountsWide3[,CallAssignements]))
write.table(CallAssignements2, "clipboard", sep="\t", row.names=FALSE)



